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On  the  bifurcation  of  the  thermo-elastic  deformation  of  an  asymmetric  plate. 
George  Jefferson,  Triplicane  A.  Parthasarathy 
UES,  Inc.,  Dayton  OH  45432 


Abstract 

The  deformation  of  a  thin  flat  elastic  plate  typically  exhibits  classic  elastic  bifurcation  behavior  when 
subjected  to  a  through-thickness-graded  thermal,  or  other  inelastic,  strain.  For  small  strain  there  is  a 
unique  solution  to  the  governing  equations,  while  for  larger  strain  multiple  equilibrium  solutions  exist. 
Physically,  the  plate  will  initially  deform  to  a  particular  shape  or  mode  and  when  a  critical  strain  is 
reached  the  mode  of  deformation  will  snap  to  a  secondary  form.  This  problem  has  been  examined 
numerous  times  in  the  literature.  However,  in  the  special  case  of  square  symmetry,  in  both  plate  shape 
and  materials  properties,  additional  equilibrium  deformation  modes  become  available  which  are  not 
predicted  by  the  existing  calculations.  In  this  paper,  analytic  solutions  are  presented  for  these  additional 
bifurcation  modes,  and  the  stability  of  these  solutions  are  evaluated  numerically.  It  is  shown  that  for 
certain  plate  materials  and  geometry  the  special-case  modes  occur  at  a  significantly  lower  strain  than 
predicted  by  existing  analyses.  The  predictions  are  verified  by  comparison  with  a  finite  element 
simulation. 

Introduction 

When  a  nominally  flat  elastic  plate  is  subject  to  a  through- thickness-graded  intrinsic  stressing  the  plate 
will  develop  a  curvature  or  warp.  Examples  of  such  intrinsic  stressing  can  be  the  result  of  a  grading  of 
the  thermal  expansion  of  the  material  combined  with  a  uniform  temperature  change,  grading  of  the 
temperature  itself  over  a  uniform  material,  residual  stressing  resulting  from  deposition  or  other  materials 
processing  effects,  hygroscopic  effects  and  so  on.  As  the  motivation  for  this  work  was  an  experimental 
observation  of  a  thermal-expansion-graded  plate  we  will  use  the  term  thermal  strain  interchangeably 
with  inelastic  or  intrinsic  strain. 

Typically,  for  small  imposed  thermal  strain  there  is  a  unique  deformation  mode,  but  for  sufficiently 
large  strain  multiple  equilibrium  solutions  become  available.  At  the  critical  point  where  multiple 
solutions  become  available  the  small  deformation  'primary'  solution  looses  stability  and  the  deformed 
shape  of  the  plate  will  transition  to  one  of  the  secondary  solutions,  i.e.  a  stable-symmetric  bifurcation  is 
observed  (Thompson  and  Hunt,  1973).  For  example,  in  the  case  of  in-plane  isotropy,  in  both  the  elastic 
response  and  intrinsic  driving  stress,  the  initial  deformation  shape  is  bowl-like,  i.e.  a  locally  spherical, 
shape  develops,  irrespective  of  the  shape  and  dimensions  of  the  plate.  For  larger  stressing  however  the 
deflection  is  approximately  ellipsoidal  and  the  plate  has  a  cylindrical  appearance  when  the  ellipse  aspect 
ratio  is  large.  In  this  paper  for  we  will  refer  to  this  as  a  cylindrical  deformation.  The  transition  point 
from  spherical  to  cylindrical  is  dependant  on  the  dimensions  of  the  plate.  The  orientation  of  the  axis  of 
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the  cylindrical  deformation  is  determined  by  the  shape  of  the  plate,  as  well  as  the  constituent  response, 
i.e.  for  isotropic  materials,  the  orientation  is  arbitrary  for  a  circular  plate  and  oriented  parallel  to  one 
edge  for  a  rectangular  plate.  Figure  1  shows  the  various  deformation  shapes  schematically. 

This  problem  has  been  extensively  examined  in  the  literature.  Hyer  (1982)  first  determined  the  solution 
for  a  rectangular  cross-ply  laminate.  A  cross-ply  laminate  is  fabricated  by  stacking  uniaxial  fiber 
reinforced  composite  layers  in  orthogonal  directions  in  the  upper  and  lower  halves  of  the  plate.  Because 
of  the  substantial  anisotropy  in  the  thermal  expansion  of  the  laminae,  thermal  expansion  induced 
deformation  is  initially  saddle-shaped,  and  exhibits  a  bifurcation  into  a  cylindrical  form.  Notably  the 
cylindrical  secondary  solution  here  is  the  same  general  shape  as  in  the  isotropic  case.  Masters  and 
Salamon  (1993)  modified  the  Hyer  result  by  introducing  additional  degrees  of  freedom  into  the  assumed 
deformation  form,  resulting  in  an  improved  prediction  of  the  bifurcation  point  without  significantly 
increasing  the  complexity  of  the  calculation.  Freund  (2000)  extended  these  results  to  the  problem  of  a 
thin  film  on  a  circular  plate. 

Recently,  thermally  induced  plate  warping  was  again  observed  experimentally  and  compared  with  the 
known  solutions  (Parthasarathy,  Keller,  et  al ,  2004).  In  this  case,  square  plates  of  ceramic-fiber  woven 
cloth  embedded  in  a  ceramic  matrix  were  fabricated.  The  cloth  contained  equal  number  of  fibers  in  each 
direction  and  the  composites  were  fabricated  with  a  '0/90'  layup.  Fibers  with  different  thermal  expansion 
and  elastic  response  were  used  in  the  upper  and  lower  half  of  the  plate.  On  cooling  from  processing 
temperatures  the  plates  warped  to  a  cylindrical  shape  for  'large'  plates  (Figure  2)  and  a  spherical  shape 
for  'smaller'  samples.  Note  in  this  case  the  elastic  response  is  square-symmetric  (within  each  layer  the 
elastic  moduli  are  equal  in  each  direction  parallel  to  the  plate  edges).  It  is  not,  however,  isotropic  as  the 
modulus  in  the  45  degree  orientations  is  lower  than  in  the  0/90  directions.  The  reduction  in  shear 
modulus  in  this  example  is  small  so  that  isotropy  was  assumed  in  comparing  the  experiment  with 
analytic  results.  The  thermal  straining  is  symmetric  and  uniform  in  each  layer  and  the  spherical  small 
strain  result  that  was  observed  is  in  fact  what  is  predicted  analytically. 

Interestingly,  however,  while  the  observed  cylindrical  large-strain  mode  appears  to  be  in  agreement 
with  the  analyses,  the  cylindrical  bifurcation  mode  differs  subtly  from  the  predictions.  Note  in  Figure  2b 
that,  while  the  plate  curvature  is  cylindrical,  the  cylinder  axis  is  aligned  with  a  plate  diagonal  rather  than 
a  plate  edge.  This  deformation  mode  clearly  exists  in  this  case  because,  owing  to  the  square  symmetry  of 
the  plate,  the  plate  diagonals  also  represent  symmetry  planes  for  the  deformation. 

The  present  analysis  is  a  modification  of  the  Masters  and  Salamon  analysis  procedure  that  reveals  the 
'new'  bifurcation  mode  and  shows  clearly  why  this  is  observed  only  for  an  orthotropic  plate.  The  paper  is 
organized  as  follows:  first  we  review  in  some  detail  the  existing  models,  then  we  present  the  new 
analytic  solution.  Because  the  analytic  results  only  predict  the  existence  of  multiple  equilibrium 
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solutions,  in  the  following  section  the  stability  of  the  various  solution  is  examined  numerically.  Next  a 
finite  element  simulation  is  presented  for  comparison  and  verification,  and  finally  the  extension  of  the 
analytic  method  to  high  order  polynomial  forms  is  investigated. 

Parallel  mode  bifurcation 

The  bifurcation  analysis  proceeds  by  development  of  approximate  displacement  functions,  calculation  of 
the  associated  elastic  strain  energy  and  application  of  the  minimum  potential  energy  theorem  to 
determine  if  one  or  more  equilibrium  solutions  exists.  Note  that  the  results  are  approximate  and  limited 
to  the  function  space  of  the  assumed  displacement  functions.  In  order  to  obtain  closed  form  results  the 
simplest  possible  displacement  functions  are  desired,  and  so  it  is  helpful  to  know  a  priori  the  form  of  the 
displacement  and  construct  displacement  functions  that  capture  the  expected  result.  However  the  present 
analysis  is  an  example  where  an  unexpected  bifurcation  mode  is  missed  as  a  result  of  the  selection  of 
approximating  functions 

Both  Hyer  and  Masters  and  Salamon  analyze  the  plate  problem  by  assuming  an  out-of-plane 
displacement  wp[x,y)  of  the  form, 

^p(^y)=^-(ax2  +by2)  (!) 

in 

where  h  is  the  plate  thickness,  x  andy  are  the  in-plane  coordinates  in  an  axis  system  centered  at  the  plate 
centroid  and  aligned  with  the  plate  edges  (or  with  the  principal  axes  of  the  material  orthotropy  when 
considering  circular  plates).  The  dimensionless  constants  a  and  b  are  to  be  determined.  By  construction, 
a/h  and  b/h  are  the  principal  curvatures  at  the  center  of  the  plate.  Clearly,  this  form  admits  the  spherical 
mode  with,  a  =  b  the  saddle  mode  with  a  =  -b ,  and  cylindrical  modes  when  one  of  a  or  b  is  zero.  As 
noted,  the  "cylindrical"  mode  is  actually  more  closely  approximated  as  an  ellipse  with  a  *  b ,  and  one  of 
the  principle  curvatures  tending  to  zero  for  large  strain.  It  is  apparent  that  the  diagonally-oriented- 
cylinder  mode  was  not  anticipated  and  hence  is  precluded  from  these  solutions. 

Using  symmetry  arguments,  Masters  and  Salamon  propose  general  in-plane  displacements  of  the  form, 

up  (x,y)  =  df  x  +  dp  x3  +  dpxy2  ... 

(x,y)  =  dpy  +  dpy3  +  dpyx2  (  ’ 

where  up  and  vp  are  the  x-  and  y-displacements  of  the  plate  mid-plane,  and  the  d f  represent  six 
additional  constants  to  be  determined.  Note  that  a  and  b  have  been  eliminated  from  the  in-plane 
displacement  form  used  by  Masters  and  Salamon  without  loss  of  generality. 

The  strains  throughout  the  plate  are  then  calculated  using  Kirchoff  plate  theory  (Vinson  and  Chou, 
1975), 
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where  z  is  the  thickness  coordinate.  Hyer's  calculation  is  a  special  case  where  yxy  =  0  is  prescribed  on 
the  mid-plane  (z  =  0)  by  setting  df  =  d%  =  -ab/4/r .  Note  however  that  this  (ad  hoc)  constraint 
equation  ultimately  results  in  only  a  minimal  algebraic  simplification  when  compared  to  the  full  (i.e. 
Masters  and  Salamon's  "case  3")  solution.  Indeed  closed  form  solutions  may  be  obtained  with  even  more 
degrees  of  freedom  introduced  in  the  in-plane  displacements,  as  will  be  discussed  in  a  later  section. 

The  material  constitutive  response  is  written  in  'reduced'  composite  laminate  notation  (Vinson  and  Chou, 
1975), 
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where  the  material  constants  Qtj  and  the  thermal  strains  ,  £y  ,are  variable  through  the  thickness  but 
constant  over  the  plate  area.  The  strain  energy  is  then, 


'*'=  {  jjQiik- 4  J  +\o22 (s'y  - £y  J  +|q,2 (ex  - 4 \ey  - 4y )+  2Q66y2xydAdz  (5) 

-h/2AZ  1  l 

where  the  inner  integral  is  taken  over  the  plate  area  A.  Because  the  Qtj  are  functions  of  z  only  this  can  be 
integrated  explicitly,  in  terms  of  the  plate  stiffness  matrices  (A,B,D)  which  are  given  in  the  appendix. 
The  theory  of  minimum  potential  energy  requires  that,  for  an  equilibrium  solution,  &¥  =  0  (Sokolnikoff, 
1956),  hence, 


~da 


d¥ 

~ddf 


=  0 


^  =  0 

(6a) 

db 

i  =  1,6 

(6b) 

Expression  (6b)  is  a  system  of  six  equations  which  are  linear  in  the  six  df ,  and  quadratic  in  a  and  b, 
hence  the  df  may  be  directly  expressed  in  terms  of  a  and  b.  Upon  substitution,  eq  (6a)  is  then  a  pair  of 
nonlinear  expressions  in  a  and  b  only, 
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G,b  +  hG2  ,  H,a  +  hH2 
a  = - - - £ —  b  = - - - £ — 

G3(b/hf+G4  H3  {a  /  h)2  +  H4 

where  G, ,  H,  are  the  algebraic  functions  of  the  material  properties  and  plate  dimensions  given  in  Masters 
and  Salamon  (1993).  The  system  has  either  one  or  three  real  solutions  corresponding  to  the  primary  and 
secondary  deformation  modes. 

At  this  point  we  will  specialize  to  the  case  of  square  symmetry.  That  is,  the  material  constants  are  all 
symmetric  with  respect  to  a  1-2  transformation  ( 0n  =  Q22 ,  etc),  the  thermal  strain  is  symmetric, 
£y  =  £*x  =  and  the  geometry  is  square-symmetric.  Note  that  the  plate  need  not  be  a  square,  i.e.  a 
circular  plate  also  has  the  required  symmetry.  The  material  response  may  be  non-isotropic  with 
066  *  (fin  -0i2  )/2  •  For  this  symmetry,  Masters  and  Salamon's  "isotropic  simplification"  may  be 
utilized  directly,  noting  that  they  only  require  the  1-2  symmetry.  That  is,  their  result  is  more  general  then 
they  asserted.  For  this  case,  G,  =  H,  for  all  i  in  (7).  The  expressions  thus  may  be  combined  to  form  a 
fifth  order  polynomial  in  a, 

(i H£2+(h *  +Hl)a  +  L*4Hy\Hl(Hl-Hl)/(H*3L*4)+Hs2a+H*4a2)=  0  (8) 

and  an  identical  polynomial  in  b.  This  is  the  result  obtained  by  Masters  and  Salamon,  however  for 
convenience  we  have  defined  H*  and  H2  so  that  all  of  the  "H"  constants  have  equal  dimension  and 
depend  only  on  the  section  properties  and  the  shape  (not  size)  of  the  plate.  L *  is  the  dimensionless  size 
of  the  plate;  L*  =  LI  h  for  a  square  with  edge  length  L  and  is  likewise  L*  is  the  normalized  diameter  for 
a  circular  plate.  The  superscript,  s,  has  been  introduced  to  indicate  that  only  H2  depends  on  the  thermal 
strain.  Expressions  for  the  square- symmetry  H* ,  for  square  and  circular  plates,  are  given  in  the 
appendix. 

The  cubic  factor  in  Eq.  (8)  always  has  one  real  root  corresponding  to  the  a  =  b  primary  solution,  noting 
that  the  same  cubic  is  obtained  by  setting  a  =  b  in  Eq.  (7).  For  typical  stable  materials  and  grading, 
H*3  >  0  and  H*4  >  H*  >  0  and  it  is  readily  shown  that  there  is  exactly  one  root. 

The  quadratic  factor  has  two  real  roots  only  if 

(c-h‘J>4(h;J(h;-h;)/h;  (9) 

Where  the  equality  is  satisfied,  Eq  (8)  has  a  triple  root  at  a  =  -H2 /  2/7* .  This  is  the  loss  of  uniqueness, 

or  bifurcation  point.  The  right  side  of  (9)  is  notably  dependant  only  on  the  shape  of  the  plate  and  the 
section  properties.  The  thermal  strain  only  appears  in  H2  ,  hence  the  critical  thermal  strain  required  for 
the  loss  of  uniqueness  is  given  in  closed  form.  For  larger  L*2H2  the  quadratic  formula  yields  two  real 
solutions  which  are  the  major  and  minor  axis  of  the  elliptical  deflection.  In  the  large  strain  limit  the 
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curvatures  asymptote  to  0  and  a  =  -  H2  / H\  .  The  elastic  stability  of  the  solutions  is  discussed  in  a 
subsequent  section. 


Diagonal  mode  bifurcation 

As  noted,  the  assumed  deformation  form,  Eqn's  (1)  and  (2),  permits  only  the  primary  and  edge-parallel 
cylinder  modes.  In  order  to  examine  the  diagonal  mode  we  assume  displacement  forms  that  admit  the 
observed  shape.  Equation  (1)  is  replaced  with, 

wd{x,y)  =  ^-(a{x  +  y)2 +b{x-y)2)  (10) 

4  h 


where  the  constants,  a  and  b,  to  be  determined,  again  are  the  dimensionless  principal  curvatures  of  the 
deformed  shape.  In  this  case,  of  course,  the  principal  curvature  axes  are  aligned  with  the  plate  diagonals. 
Note  that  here  we  have  obviously  precluded  the  edge-parallel  bifurcation  modes,  however,  as  with  the 
parallel  mode  form,  a  =  b  corresponds  to  the  primary  'bowl'  shape.  This  approach  is  followed  because 
solving  the  diagonal  mode  problem  separately  leads  to  closed  form  results,  while  solving  the 
minimization  problem  for  a  more  general  displacement  form,  with  an  additional  constant  that  would 
permit  both  modes,  does  not. 

It  may  be  noticed  that  the  'diagonal'  bifurcation  solution  could  equivalently  be  obtained  as  a  'parallel' 
bifurcation  of  a  diagonally-oriented  plate,  by  transforming  the  Q  matrix  as  well  as  the  area  of  integration 
in  Eq  (5).  The  approach  adopted  here  is  no  more  cumbersome  and  will  facilitate  investigation  of  the 
stability  of  the  combined-mode  problem. 

Noting  that  we  are  here  considering  only  the  square  symmetry  case,  for  diagonal  displacement 
symmetry,  the  in-plane  displacements  must  satisfy,  u(x,y)  =  -u(-  x-y)  and  v(x,y)  =  u(y,x)  hence  the 
trial  displacement  form  used  is, 

ud  (x,y)  =  dd  x  +  d2dx3  +  d3dxy 2  +  ddAy  +  ddy 3  +  dd6yx2 
vd  (x,y)=  dd  y  +  ddy3  +ddyx2  +ddx  +  ddx3  +ddxy2 


Note  that  we  have  the  same  number  of  unknowns  as  in  the  parallel-mode  analysis.  Following  the 
identical  procedure  we  arrive  at  essentially  the  same  form  as  Eqn's  (7)  and  (8) 

(h;~x)c+h- 

i,4»,v +(»;+; 1)  b  (  } 

(h2s  +  (h[  +  H'l)a  +  Cny  \(h*4  +  a\h*4  -  Hi  +  2/t  )/(tf  3*Z,*4 )+  H£2a  +  (tf4*  +  A)a2)=0  (12b) 


Where  the  H*  are  the  same  as  in  Eq  (8),  and 
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(13) 


X  = 


{A\  i  +  Al2  ) 


2  B, 


66 


rBu  -B12'2 


2  B, 


An  A 


12 


66  J 


2  A, 


66 


+  (A]  i  a12)d{ 


66 


1- 


Ai  A2 


2D, 


66  J 


recalling  that  A,B,D  are  the  plate  stiffness  matrices  as-defined  in  the  appendix.  Clearly  then,  comparing 
(8)  and  (12b),  the  cubic  factor  is  the  same  and  is  likewise  obtained  by  setting  a=b  in  (12a)  so  that  the 
primary  solution  displacement  field  is  identical  to  that  of  the  parallel  mode  analysis.  Multiple  solutions 
exist  when  the  quadratic  part  of  12b  has  real  roots  hence, 

(l*2H£2  J  >  4 (h*4  + xf  (h*4  -  Hi  +  2A)/h;  (14) 

In  order  to  illustrate  the  physical  significance  of  A ,  assume  that  throughout  the  plate  the  in-plane  shear 
stiffness  is  given  by, 

066  =  ~'(011  —  012  )  (15) 


where  (j)  is  an  arbitrary  constant,  with  (j>  =  1  resulting  in  isotropy,  then  we  will  have,  from  Eq.  (Al) 
Am,  =  Pi 1  —  d|2  )/2 ,  B66  =  (f){B[  j  —  Bu )/2  ,  and  D66  =  l  —  Dn  )/2  so  that, 

aJ*-1)A,+4i  [(4, -An\D„ -D,2)-(B„ -g,2)J]  (16) 

Thus  A  can  be  seen  as  an  indicator  of  isotropy;  negative  where  the  plate  shear  modulus  is  lower  than  the 
equivalent  isotropic  value  and  positive  if  the  plate  has  high  in-plane-shear  stiffness.  (The  quantity  in  [  ] 
is  strictly  positive  for  stable  materials).When  the  elastic  response  is  isotropic,  X  =  0  so  that  Eq  (8)  and 
Eq(12b)  are  identical  and  predicted  bifurcation  points  are  the  same  for  both  modes.  As  will  be  discussed 
in  the  simulation  section,  however,  this  exact  coincidence  of  the  two  solutions  is  an  artifact  of  the 
approximating  functions  selected. 

For  nonzero  X  both  parallel  and  diagonal  modes  must  be  considered.  A  loss  of  uniqueness  occurs  when 
either  of  the  conditions  (9)  or  (14)  is  satisfied.  When  both  inequalities  are  satisfied  there  are  five 
equilibrium  solutions.  Comparing  Eq.  (9)  and  Eq.  (14)  it  can  easily  be  shown  that  so  long  as  H]  >  0  and 
H*4  >  Hi  >  0  then  a  negative  X  (i.e.  low  in-plane  shear  stiffness)  results  in  the  diagonal  mode  occurring 
at  a  lower  thermal  strain  and  a  positive  X  results  in  the  parallel  mode  occurring  first. 


Example  and  comparison  with  experiment 

Oxide  fibers  offer  a  convenient  approach  to  validate  the  model  developed  in  this  work.  In 
particular,  oxide  fibers  from  3M,  Inc.,  range  in  modulus  from  about  150  GPa  to  380  GPa  and  thermal 
expansion  coefficient  from  4  to  8.3  ppm/C.  Among  these  fibers,  Nextel610  and  Nextel720  fibers  are 
most  resistant  to  high  temperature  treatments.  Nextel610  has  a  modulus  of  380  GPa  and  CTe  of 
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8.3ppm/C,  while  Nextel720  has  a  modulus  of  262  GPa  and  a  CTE  of  6.5  ppm/C.  To  validate  the  model 
predictions,  composites  that  contained  Nextel610  fibers  on  one  half  of  its  thickness,  and  nextel720  fibers 
on  the  other  half  were  considered.  Porous  alumina  matrix  composites  were  prepared  using  2D  woven 
cloths  of  Nextel610  and  Nextel720  fibers.  They  were  processed  at  high  temperature  (~1200C)  to 
densify  the  matrix.  At  the  end  of  sintering  at  this  temperature  the  entire  composite  is  in  the  iso-strain 
condition.  Upon  cooling  to  room  temperature  the  difference  in  modulus  and  CTE  results  in  a  gradient  of 
CTE-mismatch  induced  strain.  This  strain  results  in  a  bending  of  the  composite  laminate. 

Two  size  of  composites  were  fabricated.  One  of  them  was  a  circular  plate  that  was  0.12  m 
(4,75")  in  diameter.  This  was  fabricated  using  a  layup  of  woven  (8HS)  cloths  of  the  two  fibers, 
Nextel610  and  Nextel720  (3M  corp.).  The  cloths  were  cut  into  127mm  (5")  diameter  and  stacked.  3 
plies  of  Nextel610  fiber  cloth  was  stacked  on  top  of  3  plies  of  Nextel  720  fiber  cloth.  The  cloths  were 
desized  at  600C,  lh  in  air,  before  infiltrating  with  a  slurry  of  alumina  powder  (AKP-53,  Sumitomo, 
Inc.).  Final  sintering  was  done  at  1200C,  5hrs  in  air.  The  composite  was  slow  cooled  with  the  furnace.. 

The  second  set  of  composites  were  made  to  prescribed  composition  of  fibers  by  COI,  Inc.  (now  ATK- 
COI).  The  composition  was  the  same  as  the  smaller  sample,  but  the  size  was  larger  and  the  shape  was 
square.  Plates  of  dimensions  152mm  x  152  mm  (~6"x6")  were  fabricated  and  sintered  at  an  undisclosed 
temperature,  but  close  to  1200C. 

The  as-fabricated  samples  showed  distortion  resulting  from  CTE  mismatch  upon  cooling  from 
processing  temperature.  The  circular  disc  that  was  120  mm  in  diameter  exhibited  a  symmetric  spherical 
distortion,  while  the  152  mm  square  disc  exhibited  an  asymmetrical  distortion,  with  the  undistorted  line 
following  the  diagonal  of  the  square. 

Consider  a  square  plate  of  thickness  h  =  \.25mm  and  variable  edge  length  L.  Assume  two  equal¬ 
thickness  layers,  with  elastic  constants  given  by, 


3/4£  0<z<h/2 

£  —h/2<  z  <0 


Q12  =  y  Qi 


12  -  K  W1 1 


 (/>(\-vpu 


(17) 


where  v  is  Poisson's  ratio,  £  is  the  uniaxial  stiffness  of  the  lower  lamina.  For  simplicity  of  example  we 


take  v  =  0. 1  thought  the  thickness  and  note  that  £  need  not  be  specified  as  it  factors  out  of  the  analysis. 
The  CMC  experiment  shown  in  Figure  2  is  approximated  by  setting  ^  =  0.4,  while  (f>  =  1  indicates 
isotropy. 
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Further,  assume  the  thermal  strains  are  similarly  graded  as 


f3/4«Ar  0  <  z  <  h/2 

[  a  AT  —h/2<z<0 


(17a) 


where  a  is  thermal  expansion  coefficient  of  the  lower  plate  half,  taken  to  be  a  =  8  T  0  6  /  °C  for  this 
example  and  AT  is  a  prescribed  temperature  change. 

Figure  3  shows  the  isotropic  {</)  =  l)  solution,  for  L  =  150 mm  which  is,  since  A,  =  0 ,  the  same  for  either 
mode.  The  asymmetric  expansion  plate,  when  heated  or  cooled  through  a  uniform  temperature  change  of 
486  °C  or  greater  is  predicted  to  deform  to  a  cylindrical  shape.  The  stability  of  the  solutions,  indicated 
by  the  line  types  on  the  figure  will  be  discussed  in  a  subsequent  section. 


Figure  4  shows  the  effect  of  reducing  the  shear  modulus  of  the  plate  by  setting  ^  =  0.4,  approximating 
the  CMC  experiment.  The  location  of  the  equivalent  isotropic  bifurcation  point  is  shown  for  reference. 
Reducing  the  shear  modulus  increases  the  primary  mode  deflection  slightly,  as  A66  appears  in  the 
constant  H]  ,  however  the  affect  on  the  bifurcation  solutions  is  far  more  significant.  The  edge-parallel 
mode  shifts  upward  by  ~  100  °C  and  the  diagonal  mode  shifts  downward  nearly  ~  300  °C ,  so  that  the 
diagonal  mode  is  the  one  that  will  appear  first. 


Figure  5  shows  the  sensitivity  of  the  critical  point  to  the  plate  orthotropy.  For  a  plate  where  the  in-plane 
shear  is  equal-to  or  greater-than  the  corresponding  isotropic  shear  stiffness  the  lower  (i.e.  design 
limiting)  bifurcation  is  always  the  parallel  one  and  further  there  is  very  little  sensitivity  to  increasing 
stiffness.  Hence  the  isotropic  assumption  is  a  good  approximation.  On  the  other  hand,  for  reduced  shear, 
as  is  common  with  "0/90"  composites,  the  critical  thermal  strain  resulting  from  the  diagonal  mode  drops 
off  sharply,  hence  the  isotropic  assumption  may  seriously  overestimate  the  critical  condition. 


Figure  6  shows  the  trend  in  the  solutions  with  the  size  of  the  plate.  The  plate  size  appears  only  in  the  left 
term  of  Eq.  (9)  and  Eq.  (14)  and,  in  this  example,  H\  is  linear  in  the  thermal  strain,  hence  the  critical 
strain  varies  simply  as  (h/L  f  .  From  this  figure  we  can  identify  the  critical  plate  size  where  instability 
occurs.  Figure  6  also  shows  the  critical  temperature  variation  with  diameter  for  circular  plates.  For  the 
same  dimension,  diameter  or  edge  length,  circular  plates  are  somewhat  more  stable. 


If  we  consider  now  the  motivating  experimental  results  we  can  show  only  an  approximate  correlation. 
The  model  correctly  predicts  a  diagonal  deformation  of  the  square  plate,  however  the  critical  size  for  a 
circular  plate  is  predicted  to  be  only  80mm  diameter.  There  are  several  reasons  for  this.  Principally  as 
will  be  discussed  in  the  simulation  section,  the  bifurcation  analysis  tends  to  be  conservative.  The  true 
critical  size  for  a  specified  thermal  strain  will  be  somewhat  larger  than  predicted.  Additionally,  the 
material  properties  are  idealized  approximations.  In  particular,  the  actual  difference  in  thermal 
expansion  between  the  plys  in  the  composite  is  likely  somewhat  smaller  than  assumed  in  the  simulation, 
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making  the  prediction  even  more  conservative.  None  the  less  the  critical  size  is  the  correct  order  of 
magnitude,  and  more  importantly  the  diagonal  form  is  correctly  predicted  for  the  square  plate. 


Stability 

The  results  obtained  thus  far  show  only  the  existence  of  multiple  solutions.  In  order  to  evaluate  the 
stability  of  the  solutions  we  must  construct  the  Hessian  matrix, 


a2T 

OPiPj 


i,j  =l,N 


(18) 


where  the  p,  represent  all  of  the  N  unknowns  for  the  respective  solution  (i.e.  p,  =  a,p2  =  b,  p3  =  d, ... ) 
When  the  Hessian  matrix  is  positive  definite  the  solution  is  a  local  minima  of  the  energy  function,  hence 
a  stable  solution.  Unfortunately,  algebraically  evaluating  the  positive  definiteness  of  the  matrix  is  not 
tractable,  so  we  proceed  by  populating  the  matrix  using  the  analytic  solutions  and  then  numerically 
evaluating  determinants  of  all  of  its  principal  minors.  Using  this  procedure  both  the  parallel  and  diagonal 
solutions,  considered  separately,  exhibit  the  expected  'exchange  of  stability'  behavior  typical  of  such 
elastic  bifurcations  (Thompson  and  Hunt,  1973).  That  is,  the  primary  mode  is  stable  up  to  the  bifurcation 
point,  beyond  which  the  primary  mode  looses  stability  and  each  of  the  complementary  asymmetric 
solutions  is  stable.  Its  important  to  note  however  that  this  stability  determination  is  valid  only  for 
deformation  restricted  to  the  assumed  function  space.  In  other  words,  it  is  possible  that  an  unstable 
configuration  is  erroneously  predicted  to  be  stable  because  we  not  have  not  used  sufficiently  general 
trial  functions. 


In  order  to  address  the  combined  stability  consider  a  displacement  form  which  is  sufficiently  general  to 
include  either  parallel  or  diagonal  bifurcation  modes, 

w*(x,y)  =  -*-(aV  +b*y2  +c*xy)  (19) 

2  h 

u*{x,y)  =  dlx  +  d*2x3  +  d3*xy2  +d7*y  +  c/*y3 +d9‘yx2 

v*(x,y)=  d*4y  +  d5*y3  +d*6yx2  +  d7x  +  d8*x3  +  d9xy2 

*  /  *  *  *  \  * 

and  let  'T  =  T(w  ,u  ,v  j  be  the  associated  potential  energy  functional.  Stationary  solutions  for  T 

are  found  by  solving  ST  j  da  =  S'?  fdb  =  ST  jdc  =  &¥  j  <9d,  =  0  to  eliminate  the  d,  as  before. 

*  *  * 

The  resulting  three  nonlinear  equations  in  a  ,b  ,c  are  not  readily  solvable,  however  it  can  be  shown 
that  each  of  the  solutions  to  the  separate-mode  problems  also  satisfies  the  stationary  condition  for  the 
combined  mode  problem,  i.e.  with  a  =  a,b  =  b  while  c  =  0  for  parallel  solution  and  c  =  a- b  for 
diagonal.  The  algebra  required  to  show  this  is  lengthy  but  straightforward.  Of  course  it  is  plausible  that 
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even  more  solutions  exist  for  the  formulation  with  three  degrees  of  freedom  in  the  out  of  plane 
deformation,  however  numerical  investigation  did  not  reveal  any  new  modes. 

For  each  of  the  known  solution,  the  Hessian  matrix  of  the  T'  function  must  be  checked  for  positive 
definiteness.  For  the  isotropic  case  the  result  is  always  the  same  as  was  found  for  the  separate  modes, 
that  is  each  of  the  asymmetric  modes  is  stable  for  any  thermal  strain  beyond  the  shared  bifurcation  point. 
Indeed  both  modes  show  the  same  strain  energy  at  any  point  along  the  curve,  and  so  there  is  nothing  to 
indicate  that  one  solution  should  be  favored.  However  as  will  be  discussed  in  the  next  section,  the  fact 
that  the  two  modes  show  precisely  the  same  bifurcation  point  is  an  artifact  of  the  chosen  displacement 
forms. 

In  the  anisotropic  case,  there  are  two  different  bifurcation  points.  As  may  be  expected  the  solution  that 
exhibits  the  lower  bifurcation  (i.e.  diagonal  in  the  \<Q  example)  exhibits  the  usual  exchange  of  stability 
and  the  bifurcation  solution  remains  stable  for  large  thermal  strains.  However,  the  mode  that  occurs 
later  (parallel)  also  looses  stability  at  the  lower  bifurcation  point  when  the  combined  displacement 
function  is  considered.  That  is,  between  the  two  bifurcation  points  the  solution  appears  to  be  stable  if  the 
existence  of  the  first  mode  is  ignored  but  is  actually  unstable.  This  condition  is  denoted  quasi-stable  on 
Figure  4.  Interestingly,  somewhat  beyond  its  bifurcation  point,  the  parallel  curvature  solution  regains 
stability,  that  is,  it  is  a  local  minima  of  the  combined  problem.  The  potential  T'  is  however  always 
lower  for  the  diagonal  mode,  making  it  the  global  minimum. 

Simulation 

In  order  to  confirm  the  results,  the  analysis  was  duplicated  using  an  ABAQUS  finite  element  model 
(ABAQUS,  2004).  Figure  7a  shows  the  primary  deformation  mode  for  the  plate,  using  a  mesh  of  1000 
shear-deformable- shell-elements  (type  S4).  As  is  typical,  the  FEM  solution  does  not  reveal  the 
instability  of  the  primary  solution  and  so  only  this  mode  is  obtained  from  a  simple  direct  analysis. 

A  common  procedure  to  investigate  bifurcation  behavior  is  to  impose  some  defect  into  the  plate  at  zero 
temperature  to  force  the  solution  to  the  asymmetric  modes  (see  e.g.  1994  Masters  and  Salamon,  1994). 
However,  in  this  case,  because  there  are  two  sets  of  solutions  it  would  not  be  possible  to  study  that 
quasi-stable  secondary  bifurcation  behavior,  thus  an  alternate  approach  is  suggested. 

The  secondary  solutions  are  obtained  by  first  applying  displacement  boundary  conditions  to  'bend'  the 
plate  into  the  expected  shape,  i.e.  to  simulate  the  parallel  mode  a  line  of  nodes  along  y=0  is  fixed  and 
the  edge  nodes  along  y=+/-L/2  are  displaced  uniformly  by  the  amount  given  by  Eq  (1)  along  with  the 
predicted  principal  curvature  at  some  temperature,  Tmax  ,  well  above  the  expected  bifurcation  point. 
Tmax  is  then  applied  uniformly  to  the  plate.  As  a  third  step  the  'seed'  displacement  conditions  are 
removed  (except  as  needed  to  eliminate  rigid  body  displacement).  The  validity  of  the  analytic  results  are 
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supported  by  the  fact  that  the  plate  remains  in  the  imposed  cylindrical  shape  when  the  forcing 
displacement  conditions  are  removed.  For  Tmax  sufficiently  large  that  both  (9)  and  (14)  are  satisfied, 
either  mode  is  shown  as  a  valid  equilibrium  solution  in  this  manner.  It  is  notable  that  the  seed 
displacement  conditions  need  not  be  exact.  Indeed  by  applying  approximate  displacements  and 
observing  that  the  numerical  solution  takes  on  the  predicted  form  we  have  an  indication  (though  not 
proof)  of  the  stability  of  the  configuration.  As  predicted  analytically,  both  shapes  are  sufficiently  stable 
for  sufficiently  high  temperatures  that  any  reasonable  approximation  will  drive  the  numerical  solution 
into  those  forms.  Figures  7b  and  7c  show  the  simulated  parallel  and  diagonal  cylinder  configurations. 

From  the  initial  high  temperature,  the  temperature  is  ramped  down  in  small  increments  to  zero.  For  each 
increment,  the  principle  curvatures  are  determined  by  least  squares  fits  of  the  appropriate  form,  Eq.  (1) 
or  Eq  (10),  to  the  nodal  displacements  derived  from  the  finite  element  analysis. 

Figure  8  shows  the  simulated  results  using  the  same  material  parameters  as  in  the  previous  isotropic 
example,  hence  may  be  compared  directly  with  Figure  3.  The  parallel  mode  shows  reasonable  agreement 
for  the  isotropic  case,  with  the  difference  between  analytic  and  numerical  results  comparable  to  that 
found  by  Masters  and  Salamon  (1993). 


The  numerically-derived  diagonal-mode  bifurcation  point  is,  however,  higher  than  the  parallel,  thus 
contradicting  the  analytic  prediction  that  the  two  cases  should  be  identical.  Qualitative  insight  into  this 
issue  is  gleaned  by 

Note  that  the  finite  element  results  are,  of  course,  also  approximations.  However,  the  significant 
difference  in  the  methods  lies  in  the  fact  that  the  analytic  results  are  restricted  to  quadratic  out-of-plane 
approximate  displacement  while  the  finite  element  method  can  approximate  an  arbitrary  form  using  a 
much  greater  number  of  degrees  of  freedom.  By  comparing  the  numerical  solutions  with  the  least- 
squares  quadratic  form  fits  that  were  used  to  determine  their  curvatures,  we  find  that  in  for  both  modes 
the  fits  are  qualitatively  quite  good  —  that  is  the  two  modes  are  indeed  well  approximated  by  the  forms 
(1)  and  (10).  However  the  FEM  solution  reveals  deviation  from  the  simple  polynomial  displacement 
form,  especially  near  the  edges  and  comers  of  the  square  plate.  Further,  the  deviation  from  the  quadratic 
from  is  more  significant  for  the  diagonal  mode  than  for  the  parallel  mode.  On  this  bases  we  will 
conclude  that  the  FEM  predictions  are  the  most  accurate,  and  that  between  the  two  analytic  solutions  we 
expect  that  the  parallel  mode  solution  is  somewhat  more  accurate. 

We  can  conclude  then,  as  noted  previously,  that  the  'exact'  coincidence  of  the  analytic  predictions  for  an 
isotropic  plate  is  an  artificial  result  of  one  solution  using  a  marginally  more  accurate  trial  function  than 
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the  other.  This  conclusion  is  further  supported  by  the  fact  that  no  evidence  has  been  found  in  the 
literature  of  an  observation  of  a  diagonal  bifurcation  in  an  isotropic  plate. 

The  'new'  solution  is,  nonetheless,  quite  valuable  when  any  anisotropy  is  taken  into  account.  Figure  9 
shows  the  finite  element  simulation  for  a  low  in-plane  shear  stiffness  plate  (<|)=0.4).  The  numerical 
simulation  shows  the  same  shifting  of  the  two  solutions'  bifurcation  points  in  opposite  directions  as  was 
predicted  by  the  analysis. 


Higher  order  improved  accuracy  solutions. 

In  an  effort  to  close  the  difference  between  the  finite  element  and  analytic  results,  improved  accuracy 
solutions  are  obtained  by  using  higher  order  polynomials  in  the  trial  functions.  Such  higher  order 
solutions  were  obtained  for  the  parallel  mode  bifurcation  by  Salamon  and  Masters  (1994)  using  6th 
order  polynomials  in  the  strain,  or  equivalently  7th  order  displacement  functions.  Unfortunately,  as  was 
noted  in  the  stability  analysis  section,  introducing  more  than  two  degrees  of  freedom  in  the  out  of  plane 
displacement  function  results  in  an  unwieldy  system  or  nonlinear  equations.  The  in-plane  displacements 
u(x,y)  and  v(x,y)  may  in,  however  be  expanded  in  arbitrary  polynomials  in  x  and  y.  For  any  polynomial 
form,  all  of  the  unknown  coefficients  ultimately  appear  linearly  in  the  variational  statement,  hence  they 
may  always  be  (in  principle)  algebraically  eliminated  to  produce  a  pair  of  nonlinear  equations  in  the 
principal  curvatures  (a  and  b ).  Here  we  examine  the  convergence  of  the  solutions  for  increasing  order  of 
in-plane  displacement  functions. 


For  the  parallel  mode  deformation,  using  symmetry  arguments,  we  require  only  odd  powers  of  x  and 
even  powers  of  y  in  the  x-  displacement  and  likewise  odd  powers  of  y  and  even  powers  of  x  in  the  y- 
displacement.  The  following  construction  results  in  polynomial  expressions  containing  all  such  terms  up 
to  a  power  of  2 N  + 1 . 


£  £ W2,*VJ 


7=0, N  i=0,N-j 


’( x>y)=  X 


2/+1  2  j 


7=0,  N  i=0,N-j 


(20a) 


where  Ip  is  an  index  function  constructed  to  assign  unique  indices  to  each  of  the  coefficients,  / , 

Ip{i,j)  =  j{lN  +  l>-j)+2i  +  l  (20b) 

Expansion  of  Eq.  (20a)  results  in  (N+  l)(N+2)  unknowns.  Note  that  N=1  produces  trial  displacement 
forms  equivalent  to  Eq(2)  and  N=3  results  in  Masters  and  Salamon's  "expanded"  order  formulation. 

For  the  diagonal  mode,  again  using  symmetry  arguments,  we  obtain  the  expressions, 
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(20c) 


ud(x,y)=  X 

j=l,N+l  i=l,2j 

vd(x,y)=  X  Eg/d(/,y)yMx2y-/ 

y=i,A/+i  #=i,2  j 

id(i,j)=j(i-i)+i 

In  this  case,  each  displacement  function  has  (N+l)(N+2)  terms,  however  the  coefficients,  g,  are  shared 
so  that  the  same  number  of  undetermined  coefficients  results.  Likewise  N=1  results  in  the  cubic  form 
(11). 

For  both  parallel  and  diagonal  modes,  combining  these  high  order  expansions  with  the  appropriate  out- 
of-plane  expression,  (1)  or  (10),  results  in  a  system  of  nonlinear  equations  that  may  be  (in  principle) 
algebraically  reduced  to  the  form  (8).  Improved  accuracy  solutions  are  thus  obtained.  The  resulting 
expressions  for  the  "H"  constants  become  increasingly  unwieldy  with  large  N,  however.  In  fact  for  N>3 
it  becomes  impractical  to  solve  the  linear  system  equivalent  to  Eq  (6b)  even  with  a  computer  symbolic 
algebra  package.  Hence  for  large  N  we  have  used  numerical  methods  to  generate  results  for  a  specific 
example. 

It  should  be  mentioned  that  imposing  symmetry  on  the  displacement  approximation  serves  only  to 
simplify  the  algebra  and  eliminate  rigid  body  displacement.  As  there  are  no  displacement  boundary 
conditions,  the  minimum  energy  principle  requires  only  continuous,  differentiable  displacement 
approximating  functions.  Any  asymmetric  terms,  if  included,  will  naturally  cancel  from  the  analysis. 
Further,  there  is  no  reason,  other  than  algebra  reduction,  to  impose  edge  conditions  on  the  displacement. 
For  example  imposing  zero  shear  strain  on  a  stress  free  edge  is  not  necessary,  although  such 
approximations  may  lead  to  algebraic  formulations  that  can  reasonably  be  solved  by  hand  (Hyer,1982). 
In  the  present  analysis  no  such  extra  approximations  are  made,  however  our  approach  is  essentially 
dependant  on  the  use  of  computational  symbolic  algebra  routines,  and  even,  as  noted,  numeric  methods 
for  large  N. 

Figure  10  shows  that  for  increasing  polynomial  order  the  bifurcation  predictions  converge  to  fixed 
values.  It  is  readily  apparent  that  the  seventh  order  result  is  sufficient  and  that  negligible  gains  are  made 
with  larger  N.  It  is  also  apparent  that  the  analysis  results  do  not  converge  exactly  to  the  finite  element 
derived  results.  This  is  readily  seen  as  a  result  of  the  fact  that  we  are  not  adding  additional  degrees  of 
freedom  to  the  out-of-plane  trial  functions.  From  this  we  can  conclude  that  it  is  not  worthwhile  to 
expand  the  in-plane  trial  functions  to  much  higher  order  than  the  out-of-plane  trial  functions. 

Closure 

It  has  been  shown,  for  the  special  case  of  thin  plates  with  square-symmetry,  an  additional  bifurcation 
mode  becomes  available.  Most  significantly,  for  the  case  of  '0/90'  composites,  where  the  in-plane  shear 
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modulus  is  typically  relatively  low,  the  new  mode  occurs  at  much  lower  mismatch  strain  than  predicted 
by  prior  analysis.  The  analytic  method  presented  here  provides  a  conservative  estimate  of  the  critical 
condition  as  well  as  a  procedure  for  demonstrating  the  stability  of  the  various  solutions. 
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Appendix  1 

The  components  of  the  laminate  stiffness  matrix  are  defined  as  (Vinson,  1986), 
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(Al) 


hi  2  hi  2  hi  2 

4=  |  £?4-  By=  \QyZdZ  Dy  =  j  QyZ2dZ 

-A/2  -A/2  -A/2 


The  forces  and  moments  (per  unit  width)  in  the  plate  associated  with  the  thermal  strain  for  the  special 
case  of  square  symmetry  are, 

A/2  A/2 

N‘  =  \{Qn+Qnydz  M‘=  f {Qn+Qn)slzdz  (Ala) 

-A/2  -A/2 


Square  symmetry  H  constants.  Non-superscripted  quantities  are  Masters  and  Salamon's  "isotropic  case 
3"  constants.  The  corresponding  superscripted  quantities  are  normalized  to  have  equal  dimension  (stress 
cubed).  Three  of  the  constants  are  dependant  only  on  the  section  properties. 

n't  =Hjhs  =(i/A!Xa2(4!, +AM  +%)) 

Hi  =Hjh*=  (1/A)4  (4 ,  -  At  X(4 ,  + Ai  W  -  ( B , ,  +  Ba  )N'  ]  (A2) 

Hi  =  H,  /  As  =  Hi  +  (l  /  A5  )(J„  +  4,  )[(D, ,  -  D12  )(4,  -  4i)-  («„  -  B,  J  J 

The  fourth  constant  depends  on  the  shape  of  the  plate. 

For  a  square  plate,  edge  length  L, 


H*3=H3/(h7L*4) 


_U21  4)  f_ 4  4  +WAn] 


144h 


a66 


L*  =  LI  h 

For  a  circular  plate  of  radius  R 

r.  U/-4A4-4 


H,  =  , 

3  h33&4 


.-1 


+  —  2A} 
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66 


L*  =2R/h 

For  a  circular  ring  of  outer  radius  R  and  inner  radius  R0  =  g  R 

U2, -tiJh+d+d) 


H, 


L  =2  R/h 


384 h3 


AZ  _  A Z 

-Ai - a  +  64!  -242  +12(4! 
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*66 


{  4  ) 

2  ^ 

\l-yj 

J 

-1 


(A2a) 


(A3) 


(A4) 


Appendix  2  Uniform  elasticity  case. 

For  completeness,  and  to  illustrate  the  calculation,  consider  the  problem  posed  by  Freund  (2000)  of  a 
two-layer  circular  plate  with  uniform  elastic  response  throughout.  That  is  for  all  z, 
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Ql2  =V  Q„ 


Qn 


(A5) 


where  vis  Poisson's  ratio.  The  constant  Qn  may  be  written  in  terms  of  Poisson's  ratio  and  Young's 
modulus,  E,  however  it  factors  out  of  all  of  the  relevant  expressions.  A  constant  mismatch  strain  is 
specified  in  one  layer  relative  to  the  other 


J  -  J  -  A  +  Sm 

-h/2<z<(hs-hf)/: 

s'-s’-\ 

(hs-hf)/2<z<h/2 

(A6) 


where  the  plate  thickness  h  is  divided  into  'substrate'  and  'film'  regions  of  thickness  hs ,  hf  respectively, 
£s  is  the  thermal  strain  in  the  substrate  andsm  is  an  additional,  or  mismatch,  strain  in  the  film. 

The  laminate  stiffness  and  driving  moments  are, 


4=6,*  «  "  a,=  e#*712 

,  ,  s£,hfh^ 

M'=(Qn+Q,  2)-^- 


(A7) 


The  H  coefficients  and  A  are  then, 

H*=Q3nv{ 1-v2)/12 
Hs2  =  Q3n(l-v\l  +  v)2  £mhfhs 


K=Q3  n(l-v2)/12 

fa  -v2) 


2E 


H3=Qn 


768(4  +  (1  -  <p)(v  -  3)) 
A  =  fa^-Q\x(\  +  v\\-v)2 


(A8) 


so  that,  after  substitution  into  ()  the  mismatch  strains  required  for  parallel  or  diagonal  bifurcation  are, 
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r  a  V 

r  a2  ] 

+ 

1 

1 

3(1  +  vf2 

u. 

UaJ 

V  * 

(A9) 

for  (j>  =  1  both  are  equal,  and  in  the  thin  film  limit,  h  -A  hs  so  that  the  critical  strain  is, 


Sp=£d  = 


3(1  +  vj 


3/2 


V2 


(A  10) 
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which  agrees  exactly  with  Freund's  prediction  found  by  combining  Eq  (10)  and  Eq  (5)  in  that  paper. 
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spherical 


saddle 


cylindrical  (parallel)  cylindrical  (diagonal) 

Figure  1.  a)  Schematic  of  deformation  shapes.  The  spherical  shape  is  observed  for  small  strain  for  the 
square-symmetry  case.  Saddle  modes  are  observed  for  orthogonal  cross-ply  laminates.  The  edge -parallel 
and  diagonal  cylindrical  shapes  are  observed  for  large  strain  in  the  square  symmetry  case.. 
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Figure  2  Mixed  fiber  CMC  places  exhibiting  different  bifurcation  modes  a)  120  mm  circular  plate 
shows  symmetric  deformation  b)  150  mm  square  plate  shows  diagonal  deformation. 
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Figure  3.  Isotropic  plate  deformation.  A  1.25  mm  thick  by  150mm  square  plate  with  properties  given  by 

Eq  (17)  and  (17a)  with  </>  =  1. 
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Figure  4.  Reduced  shear  plate  deformation. 
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Critical  A  Temperature,  °C  Critical  A  Temperature, 


Plate  Dimension,  edge  length  L  or  diameter  2R 


Figure  6  critical  temp  vs.  plate  size. 
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Figure  7.  Finite  element  simulation,  a)  primary  solution  b)  edge-parallel  mode  c)  diagonal  mode 
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Figure  9.  Finite  element  simulation  results  for  (|)=0.4. 
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Figure  10  Convergence  of  high-order  solutions. 
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